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Abstract 

We present SUSYXATTICE - a C++ program that can be used to simulate certain classes of supersymmetric Yang-Mills (SYM) 
theories, including the well known N = 4 SYM in four dimensions, on a flat Euclidean space-time lattice. Discretization of 
SYM theories is an old problem in lattice field theory. It has resisted solution until recently when new ideas drawn from orbifold 
constructions and topological field theories have been brought to bear on the question. The result has been the creation of a new 
class of lattice gauge theories in which the lattice action is invariant under one or more supersymmetries. The resultant theories 
are local, free of doubters and also possess exact gauge-invariance. In principle they form the basis for a truly non-perturbative 
definition of the continuum SYM theories. In the continuum limit they reproduce versions of the SYM theories formulated in terms 
of twisted fields, which on a flat space-time is just a change of the field variables. In this paper, we briefly review these ideas and 
then go on to provide the details of the C++ code. We sketch the design of the code, with particular emphasis being placed on SYM 
theories with N = (2, 2) in two dimensions and N = 4m three and four dimensions, making one-to-one comparisons between the 
essential components of the SYM theories and their corresponding counterparts appearing in the simulation code. The code may 
be used to compute several quantities associated with the SYM theories such as the Polyakov loop, mean energy, and the width of 
the scalar eigenvalue distributions. 
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1. Introduction 

The problem of formulating supersymmetric theories on lattices has a long history going back to the earliest 
days of lattice gauge theory. However, after initial efforts failed to produce useful supersymmetric lattice actions 
the topic languished for many years. Indeed a folklore developed that supersymmetry and the lattice were mutually 
incompatible. However, recently, the problem has been re-examined using new tools and ideas such as topological 
twisting |[l]|2[3ia|5l|6l|7l[ai3IIS[II][l2l[l3l[Il|l5], orbifold projection and deconstruction Qgl [iTl [111 UH |20l 
|2T1|22]|23]|M]|25], and a class of lattice models have been constructed which maintain one or more supersymmetries 
exactly at non-zero lattice spacind^ 



' Present address. 

Email addresses: smcatteraliagmail . com (Simon Catterall), anoshSlanl . gov (Anosh Joseph) 
^There exist other attempts to study various supersymmetric models on the lattice. See I26II27||^|29|[^II3 1113211^^ . 
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The availability of a supersymmetric lattice construction is clearly very exciting. For example, having a lattice 
construction of the well known N — A SYM in four-dimensions is very advantageous from the point of view of 
exploring the connection between gauge theories and string/gravitational theories. But even without this connection 
to string theory, it is clearly of great importance to be able to give a non-perturbative formulation of a supersymmetric 
theory via a lattice path integral, in the same way that one can formally define QCD as a limit of lattice QCD as the 
lattice spacing goes to zero and the box size to infinity. From a practical point of view, one can also hope that some of 
the technology of lattice field theory such as strong coupling expansions and Monte Carlo simulation can be brought 
to bear on such supersymmetric theories. 

In this paper, we will briefly describe the key ingredients that go into the lattice constructions of a variety of 
SYM theories and then provide the details of C++ code that can be used to simulate these theories. In Section|2]we 
provide the general method of twisting the supersymmetries of certain classes of SYM theories to provide twisted 
SYM theories that are compatible with discretization on the lattice. We start with the twisted N - 2 SYM in two 
dimensions as a warm up example and after writing down the discretization of this theory we go on to describe the 
twisted versions of N - A- SYM in three dimensions and N - A SYM in four dimensions. In Section|3]we describe 
the algorithms used to simulate these resultant lattice theories: rational hybrid Monte Carlo (RHMC) algorithm to 
compute the fermion determinant with fractional power, leapfrog algorithm to evolve the system of equations and 
then a Metropolis test to accept or reject the configurations. In Section|4]we provide the overall structure of the C++ 
code and describe how the code advances by generating new configurations using RHMC algorithm, saves the field 
configurations after some number of Monte Carlo sweeps and measures the observables in the theory. We provide 
some simulation results in Section|5] specific to the N - 2 SYM in two-dimensions. We compute the eigenvalues of 
the scalars of the theory, study the Pfaffian phases and the presence of fermionic sign problem in the theory and also 
investigate the restoration of supersymmetry by checking if the scalar supersymmetry has indeed been implemented 
correctly in our C++ code. We provide some conclusions and outlook in Section[6] We also provide three appendices: 
[Appendix A| details the installation of the program. Appendix B| lists the files included in SUSYXATTICE library 
with a brief description of their purpose and in Appendix C we provide a sample file with input parameters. 

We hope that this work will motivate elementary particle physicists as well as high energy computational physicists 
to pursue numerical studies of supersymmetric lattice theories in particular, the AT = 4 Yang-Mills in four dimensions. 



2. The method of topological twisting in SYM theories 

First, let us explain why discretization of supersymmetric theories resisted solution for so long. The central 
problem is that naive discretizations of continuum supersymmetric theories break supersymmetry completely and 
radiative efifects lead to a profusion of relevant supersymmetry breaking counterterms in the renormalized lattice 
action. The coefficients to these counterterms must then be carefully fine tuned as the lattice spacing is sent to zero in 
order to arrive at a supersymmetric theory in the continuum limit. In most cases this is both unnatural and practically 
impossible - particularly if the theory contains scalar fields. Of course, one might have expected problems - the 
supersymmetry algebra is an extension of the Poincare algebra, which is explicitly broken on the lattice. Specifically, 
there are no infinitesimal translation generators on a discrete space-time so that the algebra {Q,Q] - jaPa, where a is 
the space-time index, is already broken at the classical level. Equivalently, it is a straightforward exercise to show that 
a naive supersymmetry variation of a naively discretized supersymmetric theory fails to yield zero as a consequence 
of the failure of the Leibniz rule when applied to lattice diff'erence operators. In the last five years or so this problem 
has been revisited using new theoretical tools and ideas and a set of lattice models have been constructed which retain 
exactly some of the continuum supersymmetry at non-zero lattice spacing. The basic idea is to maintain a particular 
sub-algebra of the full supersymmetry algebra in the lattice theory. The hope is that this exact symmetry will constrain 
the eff'ective lattice action and protect the theory from dangerous supersymmetry violating counterterms. 

Two approaches have been pursued to produce such supersymmetric actions: one based on ideas drawn from the 
field of topological field theory JUiHQ and another pioneered by David B. Kaplan. Mithat Unsal and collaborators 
using ideas of orbifolding and deconstruction lfT8l[T9ll20l . Remarkably, these two seemingly independent approaches 
lead to the same lattice theories - see lfT2l [211 l22l [341 and the recent reviews ifTSl |35l |36l . This convergence of 
two seemingly completely different approaches to the problem leads one to suspect that the final lattice theories may 
represent essentially unique solutions to the simultaneous requirements of locality, gauge-invariance and at least one 



2 



Simon Cattemll and Anosh Joseph / Computer Physics Communications 00 (2012) 1 ^25\ 



3 



exact supersymmetxy. In this paper, we will use the language of topological twisting to discuss these supersymmetric 
lattice constructions, but the reader should remember that the orbifold methods lead to the same lattice theories. 

2.1. Twisting the supersymmetries in d dimensions 

The basic idea of twisting goes back to Witten in his seminal paper on topological field theory f37] but actually 
had been anticipated in earlier work on staggered fermions |38 1. In our context the idea is decompose the fields of the 
theory in terms of representations not of the original (Euclidean) rotational symmetry S Omtid) but a twisted rotational 
symmetry SO{d)' , which is the diagonal subgroup of this symmetry and an SO^id) subgroup of the R-symmetry of 
the theory, 

S 0{d)' = diag(5 OLorentz(fl') X S 0^{d)) . (1) 

To be explicit, consider the case where the total number of supersymmetries is Q = 2'^. In this case we can treat 
the supercharges of the twisted theory as a I'^l^ x 2''^^ matrix q. This matrix can be expanded on products of gamma 
matrices: 

q^Ql + Qaja + Qabjajb + ■■■ (2) 

The 2"' antisymmetric tensor components that arise in this basis are the twisted supercharges and satisfy a correspond- 
ing supersymmetry algebra following from the original algebra 

= 0, (3) 

{Q,Qa} = Pa, (4) 

! (5) 

The presence of the nilpotent scalar supercharge Q is most important: it is the algebra of this charge that we can hope 
to translate to the lattice. The second piece of the algebra expresses the fact that the momentum is the ^-variation of 
something, which makes plausible the statement that the energy-momentum tensor and hence the entire action can be 
written in (3-exact forn|^ Notice that an action written in such a Q-exact form is trivially invariant under the scalar 
supersymmetry, provided the latter remains nilpotent under discretization. 

The rewriting of the supercharges in terms of twisted variables can be repeated for the fermions of the theory and 
yields a set of antisymmetric tensors (77, ij/a^Xab, ■ ■ ■), which for the case of Q-l"^ matches the number of components 
of a real Kahler-Dirac field. This repackaging of the fermions of the theory into a Kahler-Dirac field is at the heart of 
how the discrete theory avoids fermion doubling as was shown by Becher, Joos and Rabin in the early days of lattice 
gauge theory 139114011 . 

It is important to recognize that the transformation to twisted variables corresponds to a simple change of variables 
in flat space - one more suitable to discretization. A true topological field theory only results when the scalar charge 
is treated as a true BRST charge and attention is restricted to states annihilated by this charge. In the language of 
the supersymmetric parent theory such a restriction corresponds to a projection to the vacua of the theory. It is not 
employed in the lattice constructions we discuss in this paper. 



2.2. A warm up example: Twisted N — 2 SYM in two dimensions 

We look at the twisted N - 2 SYM in two dimensions as a warm up example. This theory satisfies our require- 
ments for supersymmetric latticization: its R-symmetry possesses an S0{2) subgroup corresponding to rotations of 
the two degenerate Majorana fermions into each other. Explicitly the theory can be written in twisted form as 



Xal,rab+T][D„,Da]-\lld\ . (6) 



The degrees of freedom are just the twisted fermions (//, i//a,Xab) previously described and a complex gauge field 
The latter is built from the usual gauge field and the two scalars present in the untwisted theory - Ag + iB„ with 
corresponding complexified field strength Tab- 



'Actually in the case of the four-dimensional N = A there is an additional (3-closed term needed. 
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The complex covariant derivatives appearing in these expressions are defined by 

Da = da+:^a^da+Aa+ iB,, , 

Da = da+^a^da+Aa-iBa. (7) 

All fields take values in the adjoint representation of t/(A^j^ It should be noted that despite the appearance of a 
complexified connection and field strength, the theory possesses only the usual U{N) gauge-invariance corresponding 
to the real part of the gauge field. 

Notice that the original scalar fields transform as vectors under the original R-symmetry and hence become vectors 
under the twisted rotation group while the gauge fields are singlets under the R-symmetry and so remain vectors under 
twisted rotations. This structure makes possible the appearance of a complex gauge field in the twisted theory. 

The nilpotent transformations associated with Q are given explicitly by 

Q^a = Iff a 

Qllfa ^ 

QMa = 

QXab = -Tab 

Qt] ^ d 

Qd ^ (8) 
Performing the <3-variation and integrating out the auxiliary field d yields 



(9) 



To untwist the theory and verify that indeed in flat space it just corresponds to the usual theory one can do a further 
integration by parts to produce 



5 = ^ J d^xTx[-Fl, + IBaDhDbBa - [Ba,Bbf + Lf) , 



(10) 



where Fgb is the usual Yang-Mills term. It is now clear that the imaginary parts of the gauge fields Ba can now be 
given an interpretation as the scalar fields of the usual formulation. Similarly one can build spinors out of the twisted 
fermions and write the action in the manifestly Dirac form 

-D2-/B2 £>i +iBi \l (Ai 
)i -iBi D2-/B2 l\ tA2 



2.3. Discretization of the twisted N — 2, d — 2 theory 

The twisted theory described in the previous section may be discretized using the techniques developed in IIT2I 
I22II25]. (Complex) gauge fields are represented as complexified Wilson gauge links 1/0(11) = living on links 

(n, n + /7^) of a lattice, which for the moment we can think of as hypercubic. These transform in the usual way under 
U (N) lattice gauge transformations 

^ G(n)%(n)G"'"(n . (12) 

Supersymmetric invariance then implies that lAfl(n) live on the same links as Kain) and transform identically. The 
scalar fermion 77(11) is clearly most naturally associated with a site and transforms accordingly 

Tj(n) ^ G(n)7,(n)G' (n) . (13) 



The generators are taken to be anti-hermitian matrices satisfying Tr{T"T ) = -6 . 
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(n + Zii +/i2) 



(n + /2i) 



Wi(n)>i(n) 



Figure 1 . The unit cell of the two-dimensional N = 2 lattice SYM with the orientation assignments for twisted fields. 



The field ;^fai(n) is sfightly more difficult. Naturally as a 2-form it should be associated with a plaquette. In practice we 
introduce diagonal links running through the center of the plaquette and choose ;^'ai(n) to lie with opposite orientation 
along those diagonal links. This choice of orientation will be necessary to ensure gauge-invariance. Fig. 1 shows the 
unit cell of the resultant lattice theory. 

To complete the discretization we need to describe how continuum derivatives are to be replaced by difference 
operators. A natural technology for accomplishing this in the case of adjoint fields was developed many years ago and 
yields expressions for the derivative operator applied to arbitrary lattice p-forms BTI . In the case discussed here we 
need just two derivatives given by the expressions 



The lattice field strength is then given by the gauged forward difference !7ai(n) = D„ 1(b(n) and is automatically 
antisymmetric in its indices. Furthermore, it transforms Uke a lattice 2-form and yields a gauge-invariant loop on the 

lattice when contracted with Xahi^)- Similarly the covariant backward difference appearing in tlaiT^) transforms 
as a 0-form or site field and hence can be contracted with the site field 77(11) to yield a gauge-invariant expression. 

This use of forward and backward difference operators guarantees that the solutions of the theory map one-to- 
one with the solutions of the continuum theory and hence fermion doubling problems are evaded f39l|. Indeed, by 
introducing a lattice with half the lattice spacing one can map this Kahler-Dirac fermion action into the action for 
staggered fermions ll42l . Notice that, unlike the case of QCD, there is no rooting problem in this supersymmetric 
construction since the additional fermion degeneracy is already required by the continuum theory. The number of 
fermions of the twisted theory exactly fills out the components of a Kahler-Dirac field and corresponds to the taste 
degeneracy of (reduced) staggered fermions. 

2.4. Twisted N — 4 SYM in three dimensions 

The twist of AT = 4 SYM in three dimensionsjcan be most succinctly written in the form where 



2)W/i.(n) = 'Z/„(n)/i(n+/7j-/i(n)'Z/„(n+A?,), 
^rVa(n) = fainfUain) - U,(n - /7„)/„(n - /7j . 



(14) 



(15) 



S 




(16) 



'This twist of N = 4, d = 3 SYM is known as the Blau-Thompson twist 1431 . 
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Figure 2. The unit cell of the three-dimensional N = 4 lattice SYM with the orientation assignments for twisted fermionic fields. 

The fermions comprise a multiplet of p-form fields (rj, iffa,Xab, ^^oAcj^ where in three dimensions p = 0, ■ ■ ■ ,3. This 
multiple! of twisted fermions corresponds to a single Kahler-Dirac field and here possesses eight single component 
fields as expected for a theory with N - 4 supersymmetry in three dimensions. 

The imaginary parts of the complex gauge field ^„ with a = 1,2,3 appearing in this construction yield the three 
scalar fields of the conventional (untwisted) theory. Fields d and Babe are auxiliaries introduced to render the scalar 
nilpotent supersymmetry Q nilpotent off'-shell. The latter acts on the twisted fields as follows 







Q^a -- 


= 


Q^a -- 


= 


QXab -- 


-- Tab 




-- d 


Qd -- 


-- 


QBabc - 


- Babe 


QOabc -- 


-- 



Notice that this construction differs slightly from the one discussed in |44). The fermion term involving a 3-form is 
here trivially rewritten as a g-exact rather than Q-closed form. 

Doing the Q-variation, integrating out the field d and using the Bianchi identity 

eabc'^cTab = , (18) 

yields 

S = ^ ^ d\i:i{- TabTab + \Wa,m^ ' Xab^la^ b} 

-IpaDaTl-OabcDicXab]). (19) 

The terms appearing in the bosonic part of the action can then be written in the following form exposing the Ba 



'it is common in the continuum literature to replace the 2- and 3-form fields in these expressions by their Hodge duals; a second vector ij„ and 
scalar f) see, for example 1431 . 
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dependence explicitly 

TabTab = {Fab ' [Ba, Bb\)(Fab ' [Ba, Bb\) + B B ij) , 

^[©«,£)„f = -2(DM\ (20) 
where F„i, and Dg denote the usual field strength and covariant derivative depending on the real part of the connection 

2.5. Discretization of the three-dimensional N — 4 SYM theory 

The transition to the lattice from the continuum theory is similar to the case of the two-dimensional N - 2 SYM 
theory. We replace the continuum complex gauge field J{.a{x) at every point by an appropriate complexified Wilson 
link'ZYa(n) = g-^"'"), a - 1,2,3. These lattice fields are taken to be associated with unit length vectors in the coordinate 
directions a in an three-dimensional hypercubic lattice. By supersymmetry the fermion fields (/^aCn), a = 1, 2, 3 lie on 
the same oriented link as their bosonic superpartners running from n — > n + /i^. In contrast the scalar fermion 77(11) is 
associated with the site n of the lattice and the tensor fermions XabiT^), a < b - 1,2,3) with a set of diagonal face links 
running from n + /i^ + /i^ — > n. The final 3-form field 6'a/,e(ii) is then naturally placed on the body diagonal running 
from n — > n + + A'c- The unit cell and fermionic field orientations of the three-dimensional theory is given in 

Fig. 2. The construction then posits that all link fields transform as bi-fundamental fields under gauge transformations 

/7(n) ^ G(ii);7(ii)G^(n) 

^ir„(n) ^ G(ii)(/r,„(ii)G"'"(n + Ai„,) 

Xmniry) G(n-H^„, +^„);i'„,„(n)G"'"(n) (21) 

Omnqiv) ^ G(n)0;„„^(n)G^(in- //,„ -H^„ 

t/,„(ii) ^ G{n)^„{n)G\n+'fiJ 

Hjn) ^ G(n+-fl,„ynjn)GHn) 

The action of the lattice theory resembles to its continuum cousin with the one modification that the continuum 
field y[fl(x) is replaced with the Wilson link U,n{a) and the lattice field strength being defined as ;?m„(n) - ''Z/„(ii). 
Thus the supersymmetric and gauge-invariant lattice action is 



S ^ Tr(^„„(n)5^„„(n) + 77(n)£)|;V,„(n) 



n,m,n,q 

+ ^77(n)flf(n) + B„,„^(n)©^^Vm«(n)) . (22) 
The covariant difference operators appearing in these expressions are defined by ll25l 

£)i;'/,„(n) = 'W„,(n)/„,(n)-/,„(n-/7,„W«(n-A?,„) (23) 
2)W/«(n) = 'Z/,„(n)/„(n+A?„,)-/;,(ii)'W„(n+A?„) (24) 

^«Vm(n) = /„(n)^„(ii)-^™(ii-/7,„)/,„(n-/7„,) (25) 

D]^\f„gin) = f„g(n+JlJ^(„(n)-^(,„(n+'^„+'^i^^)f„,|(n) (26) 

These expressions are determined by the twin requirements that they reduce to the corresponding continuum results 
for the adjoint covariant derivative in the naive continuum limit 'ZY™ — » I^v -i- and that they transform under gauge 
transformations like the corresponding lattice link field carrying the same indices. This allows the terms in the action 
to correspond to gauge-invariant closed loops on the lattice. 

Upon following the prescription |25 1 for lattice covariant derivatives, we write down the lattice action in terms of 
the link fields 1/,„(n) and U,n{n) 



g 

n,m,n,q 



-;rm„(n)£)|„>„](n) - i](n)D„, ifr,„(n) - 0,n„,(n)D^^ x,nn](n)) . (27) 
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The bosonic part of the action is 

^ Q,m,n 
^ n,m,« 

x('W,„(ii)t/„(n + /?,„) - 1/„(n)1/„,(n + /I,,)) 

+ ^(l/,„(n)^„,(n) - ^,„(n - A?„,)%,(n - /7,„))^] , (28) 

and the fermionic part 

Sf 2 Z_i Tr |-(5m^5„r - 5,„r5„cy) 

^ n,m,n,q,r,e,f 

x[x„,„in)['Ugin)tfrrin+'fi^) - i/r,(n)'ZY,/n + ^,))] 
+77(n)(i/r„,(n)1/,„(n) - 'U„,{n-'[iJ^Jn -Jij) 

^eref{n){xre{n + 'fifyUfia) - 04 fin + /7, + 11, )x rein))] . (29) 
It is easy to see that each term in the lattice action forms a gauge-invariant loop on the lattice. 

2.6. Twisted N — 4 SYM in four dimensions 

In four dimensions the constraint that the target theory possess 16 supercharges singles out a single theory for 
which this construction can be undertaken - the N - A SYM. 

The continuum twist of AT = 4 that is the starting point of the twisted lattice construction was first written down 
by Marcus in 1995 |45| although it now plays an important role in the Geometric -Langlands program and is, hence, 
sometimes called the GL-twist ||46l . This four-dimensional twisted theory is most compactly expressed as the di- 
mensional reduction of a five-dimensional theory in which the ten (one gauge field and six scalars) bosonic fields are 
realized as the components of a complexified five-dimensional gauge field while the 16 twisted fermions naturally 
span one of the two Kahler-Dirac fields needed in five dimensions. Remarkably, the action of this theory contains a 
(3-exact piece of precisely the same form as the two dimensional theory given in ^ provided one extends the field 
labels to run now from one to five. In addition, the Marcus twist requires a new (3-closed term, which was not possible 
in the two-dimensional theory. 

•^closed — g ^ ^mnpqrXqi'-^ pXmn • (30) 

The supersymmetric invariance of this term then relies on the Bianchi identity 

^mnpqr^p^qr — • (31) 



2.7. Discretization of the four-dimensional N — 4- SYM theory 

In two and three dimensions we were able to accommodate the bosonic fields of the theory in a natural way by 
assigning them to the links of a hypercubic lattice. For the £2=16 theory this is not possible; the theory can be 
parametrized in terms of five complex gauge fields in the continuum. We are thus motivated to search for a four- 
dimensional lattice with five basis vectors fi^, a - 1, ■ ■ ■ ,5. One simple solution is to use a hypercubic lattice with an 
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additional body diagonal 

= (1,0,0,0) 

^2 = (0,1,0,0) 

A<3 = (0,0,1,0) (32) 

-fi^ = (0,0,0,1) 

= (-1,-1,-1,-1) 

The field l/g is then placed on the body diagonal link. Actually, we will indeed utilize such a hypercubic lattice 
when building the C++ data structure needed to code the resulting theory. Notice that the basis vectors sum to zero, 
consistent with the use of such a linearly dependent basis. 

However, it should also be clear that a more symmetrical choice is possible in which the five basis vectors are 
entirely equivalent and the lattice theory possesses a large point group symmetry corresponding to permutations of 
the set of basis vectors. Such a discrete structure exists in four dimensions: it is called the A*^ lattice. It is constructed 
from the set of five basis vectors 'cg pointing from the center of a four-dimensional equilateral simplex out to its 
vertices together with their inverses -e^. It is the four-dimensional analog of the two-dimensional triangular lattice. 
A specific basis for the lattice is given in the form of five lattice vectors 

/II 1 1 ^ 

( ) (33) 



V2' V6' VI2' V20^ 



= (-44. ,34, 

?3 = (0,-4,^,^) (35) 

?4 = (0,0,-4=,^) (36) 
^ V12 V20^ 

?5 = (0,0,0,-^) (37) 



The basis vectors satisfy the relations 

5 



^e,„ = 0; e,n ■ e„ = (s,„„ - ^(?,„)^(?m)v = i^^v; yu, v = 1, ■ ■ ■ ,4. (38) 

m~l ;»=1 

Notice that 5^ is a subgroup of the twisted rotation symmetry group 5(9(4)'. Furthermore, the lattice fields transform 
in reducible representations of this discrete group - for example, the vector l/o decomposes into a four component 
vector 'ZY^ and a scalar field = 2a 'i/a under and hence also under 5(9(4)' in the continuum limit. Invariance of 
the lattice theory with respect to then guarantees that the lattice theory will inherit full invariance under twisted 
rotations as the lattice spacing is sent to zero. 

Complexified Wilson gauge link variables 1/,, are then placed on these links together with their (3-superpartners 
i/To. The ten twisted fermions Xab are associated with additional diagonal links'?,, +'eh with a > b while a single 
fermion rf is placed at each lattice site. 

We can connect the basis vectors of the hypercubic lattice and the A*^ lattice through a set of linear transformations 
- see Il2ni47l . The integer-valued hypercubic lattice site vector n can be related to the physical location in space-time 
using the basis vectors 

4 4 

R-a ^(/Uy ■ n)'?v = a ^ "v^v , (39) 

v=l V=I 

where a is the lattice spacing. On using the fact that Yjm'^m - 0, we can show that a small lattice displacement of the 
form dn = //^ corresponds to a space-time translation by ide,„y. 

4 4 

dR = a^(jj.y ■ dn)7y = a ^(/?^ ■ JiJ^y = de„ . (40) 

v=l v=l 

9 
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The lattice action corresponds to a discretization of the Marcus twist on this lattice and can be represented as a set 
of traced closed bosonic and fermionic loops. It is invariant under the exact Q scalar supersymmetry, lattice gauge 
transformations and a global permutation (point group) symmetry S^, and can be proven free of fermion doubling 
problems as discussed before. The (3-exact part of the lattice action is again given by (j9| with the indices fj., v now 
labeling the five basis vectors of A*^ or equivalently its hypercubic cousin. 

Finally, it is important to note that while the true lattice in space-time is this rather complicated looking A*^ 
structure, we can represent all of the lattice fields in our theory by giving only their coordinates on the abstract 
hypercubic lattice. Indeed, since the lattice action only depends on the structure of the hypercubic lattice we will 
not need the explicit coordinates of the lattice to generate Monte Carlo configurations during the simulation. 
The explicit mapping of hypercubic coordinates to space-time coordinates in the A^ lattice is only needed when, for 
example, we want to compute spatially dependent objects such as correlation functions of fields. In this case we 
should compute distances relative to the underlying A^ lattice not its hypercubic partner. 

While the supersymmetric invariance of the (3-exact term is manifest in the lattice theory it is not immediately 



clear how to discretize the continuum <3-closed term. Remarkably, it is possible to discretize ( 30 1 in such a way that 
it is indeed exactly invariant under the twisted supersymmetry: 

1 — (-) _ 

•^closed — ~Q / . ^mnpqrXqri^ f^m /^h 

+ tlp)Dp Xmnin+flp) , (41) 

which can be seen to be supersymmetric since the lattice field strength satisfies an exact Bianchi identity BTI 

^mnpqr^ p 'J^ qr ~ ■ (42) 



3. Simulating the SYM theories: Algorithms 

Although the fields entering into these twisted descriptions appear somewhat different to the usual fields used in 
QCD the basic algorithms we use to simulate them are borrowed directly from lattice QCD; namely we integrate 
out the fermions to produce a Pfaffian which is in turn represented by the square root of a determinanj^ and can be 
simulated using the usual RHMC algorithm |48|. 

If we denote the set of twisted fermions by the field ^ - (rj, ij/fi^Xtiv) we first introduce a corresponding pseudo- 
fermion field O with action 

5pF = <!)"'" (M"'"M)" I (D, (43) 

where M - M{'l(,'l(^) is the antisymmetric twisted lattice fermion operator given, for example, in 

Integrating over the fields O will then yield (up to a possible phase) the Pfaflian of the operator M(K, 1/ ' ) as 
required. The fractional power is approximated by the partial fraction expansion 

=ao + y , (44) 

where the coefficients {c^jS,) are evaluated offline using the Remez algorithm to minimize the error in some interval 
(e. A). Typically we have used P = 15 which yields a fractional error of 0.00001 for the interval 0.0000001 1000.0, 
which conservatively covers the range we are interested in. 

Following the standard procedure, we introduce momenta {pu, Pf) conjugate to the coordinates (1/, <S>) and evolve 
the coupled system using a discrete time leapfrog algorithm according to the classical Hamiltonian 

H^SB + Spf+ pnjpv + Pq,p^ ■ (45) 



27F 



'of course this ignores a possible sign ambiguity. We return to tliis issue later when we discuss whether the phase quenched simulations we use 
suffer from a sign problem. 

**The antisymmetry is guaranteed if the fermion action is rewritten as the sum of the original terms plus their lattice transposes. 
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Notice that the bosonic actioifl 



(46) 



is real, positive semi-definite in all these theories. 
One step of the discrete time update is given by 

















- 


6^^ = 


-- (e^'P-" 




- Stp^ 




- ^'f- 








- ^'f 




- y/o 



where the forces fu and /o are given by 



6S 
"50 



(47) 

(48) 

(49) 
(50) 

(51) 
(52) 

(53) 
(54) 



and the bar denotes complex conjugation. Using the partial fraction expansion given in (44i the fermionic contribu- 
tions to these forces take the form 



/■/ ermionic _ \ ^ 



_ 6M /_ 6M 
I (51/ 



/^fermionic 
ftp 



O - ^ a/Si 



where 



tj = Ms,- 



(55) 

(56) 
(57) 



(58) 
(59) 



The latter set of sparse linear equations is solved using a multi-mass conjugate gradient (MCG) solver |49|, which 
allows for the simultaneous solution of all P systems in a single CG solve. 

At the end of one such classical trajectory the final configuration is subjected to a standard Metropolis test based 
on the Hamiltonian H. The symplectic and reversible nature of the discrete time update is then sufficient to allow 
for detailed balance to be satisfied and hence expectation values are independent of 6t. After each such trajectory the 
momenta are refreshed from the appropriate Gaussian distribution as determined by H, which renders the simulation 
ergodic. 



'From now on we interchangeably use x and n to denote the lattice site. 
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The fermionic contribution to the forces are shown below 



.fennionic ^ ^ V Q,.^t zl ^(M'^ ADF 



F xt d 



1=1 

p 



!=I 

p 



F x t / F 



/ r X I um I r x n 



^ -,. dM 



rfermionic _ PJ 



«o^(/^"'>) + c^i-gp{F^[(M'^M+/3i)-'F]) 

P 

= aoFU^^ais]. (61) 
1=1 

4. Overall structure of the C++ code 

Typically the bosons lie on the usual nearest neighbor links of a hyeprcubic lattice while the fermions occupy both 
these links and additional site, face and body diagonal links. In the case of AT = 4 in four dimensions we have to 
augment the set of boson links with one additional gauge field associated with the body diagonal link of the hypercube. 
We introduce the Lattice_Vector class to store the coordinates of the lattice sites and also the vector between sites. 
Such lattice vectors can be added or subtracted by overloading the '+' or '-' operators. These operations also respect 
the lattice boundary conditions. Associated with this class is a general function loop_over_lattice(x sites) that 
implements a loop over all lattice sites indexed by their coordinate vector; thus a simple loop looks like 

while(loop_over_lattice(x, sites)) .... 

The bosonic and pseudo fermionic fields are stored in various objects which are indexed via their lattice site 
vector and whose type corresponds directly to the tensor structure of the associated continuum field so that one finds 
C++ classes labeled Site_Field, Link_Field, Plaq_Field, Body_Field etc. in the header file utilities .h. 



(We provide the list of C++ files that goes into the code in Appendix B ) The full Kahler-Dirac field is contained 
in the class Twist_Fermion while the Gauge_Field class contains the complexified Wilson gauge link. All these 
objects are in turn built from objects of type Umatrix corresponding to complex NCDLOR x NCOLOR matrices. Simple 
arithmetric operations which overload the usual arithmetic operations are defined for manipulating these objects. 

Let us briefly describe how the code works. The general organizational structure of the code is given in Fig. 3. 
We begin with sym . cpp. It reads the input parameters such as number of sweeps (SWEEPS), number of thermalization 
steps (THERM), gap in measurements (GAP), the 't Hooft coupling (LAMBDA), etc., using functions contained in the file 
read_parcmi . cpp. It can also read in previously generated field configurations using read_in . cpp. 



The code sym. cpp performs three major tasks: 
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read_param.cpp 



readjn.cpp 



measure. cpp 




update. cpp 



write_out.cpp 




kinetic_energy.cpp 




action. cpp 




force. cpp 




evolve_fields.cpp 














MCGsolver.cpp 




fermion_forces.cpp 




MCG_solver.cpp 




force. cpp 



Figure 3. The organizational structure of the C++ code that generates and measures field configurations. 



1 . Generates new configurations using a rational hybrid Monte Carlo (RHMC) algorithm. This is accomplished 
by calling the function update (U , F) contained in update . cpp. 

2. Saves the current field configuration after some number of Monte Carlo sweeps (using the functions in write_out . cpp). 

3. Measures the observables in the theory. This is done by function calls within measure . cpp. 

Let us focus on the task of updating field configurations first. After reading the initial parameters and field 
configurations update () is called. Here we refresh the momenta p_U and p_F (using a Gaussian distribution) and 
then go to kinetic_energy . cpp to compute the kinetic energy: 

Adj(p_U)*p_U + Cjg(p_F)*p_F. 

Compare this with the first two terms in the classical Hamiltonian ( [45] l: 

P^^P^| + P^P<s> ■ 

After computing kinetic energy the boson and pseudo-fermion actions ( |45] l are computed with a call to the function 
actionO . 

The computation of the bosonic action 5^ is straightforward. In the code it is accomplished with the line 

KAPPA* [0.5*Tr(DmuUmu*DmuUmu) + 2.0*Tr(Fmunu*Adj (Fmunu))] . 

Here KAPPA is the dimensionless lattice coupling. It is defined in read.param. cpp and depends on the number of 
dimensions (D), size of the lattice (LX, LY, LZ, T) and number of colors (NCOLDR). 

The code associated with spefcific terms in the bosonic action can easily be identified with its analytic expression. 
We have 

DmuUmu(x) — > Umu(x) *Udagmu(x) -Udagmu(x-e_mu) *Umu(x-e_niu) , 
Finunu(x) — > Umu(x) *Uiiu(x+e_mu) -Unu(x) *Umu(x+e_nu) . 
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The code used to compute the fermionic part of the action is given by 

S_F = ampdeg*(Cjg(F)*F) + ZHT^ amp [n] * (Cjg(F) *sol [n] ) , 

where n runs from to DEGREE (which is equal to number of terms in the Remez approximation P), ampdeg cor- 
responds to ao, F the twisted pseudo-fermion F, Cjg(F) is F' , amp[n] is ff, and sol [n] corresponds to i, = 
(M"'"M+/?i)->F. 

Again one should compare this code with the form of the pseudo-fermion action 

Spf = aoF^'F + Sf^i a,F-^ [(M-^M+pd-'F] . 

We invoke a multi-mass conjugate gradient solver MCG_solver () given in MCG_solver . cpp to help compute the 
terms needed in the fermionic action. The MCG solver can return the solutions to (M* M + /3i)si - F for all shifts /?,-. 

Once the Hamiltonian is computed we evolve the fields along a classical trajectory. This is handled by the function 
evolve_f ields. The evolution of the fields and momenta is achieved through a leapfrog algorithm. In the first half 
step we have 

p_Umu p_Umu + . 5*DT*f _Umu 

p_F p_F + 0.5*DT*f_F 

Umu — > Umu + exp (DT*p_Umu) 
F ^ F + DT*p_F 

Immediately after computing the change in fields (Umu and F) and momenta (p_Umu and p_F), we update the forces 
by calling f orce () . The bosonic force contribution to f _Umu is given by 

f_Umu(x) — > f _Umu(x)+Umu(x)*Udagmu(x)*DmuUmu(x) 
-Umu(x) *DmuUmu(x+e_mu) *Udagmu(x) 
+2 . 0*Umu (x) *Unu (x+ejiu) * Ad j (Fmunu (x) ) 
-2 . 0*Umu(x) *Adj (Fmunu (x-ejnu) ) *Unu(x-e_nu) 

The computation of the fermionic force f _F requires first a call to the MCG solver MCG_solver(). We find 

f_F = -ampdeg*Cjg(F) - Y^HT^ amp [n] *Cjg(sol [n] ) . 

Once we have this solution an additional contribution to the gauge force coming from the pseudo-fermions is gotten 
by a call to the function f ermion_f orces () . Each fermionic term in the action yields a contribution. We provide 
a part of this code in Fig. 4. In the second half step of the leapfrog algorithm the momenta p_U and p_F are again 
updated with the new forces. These final forces are then saved for the next iteration. 

In practice, it is important to use a multi-time step integrator for this evolution |50|. In this case while the fermions 
are evolved with a time step of DT, the bosons are integrated with the time step DT/MSTEP. Provided the boson force 
is substantially larger than the fermionic contribution this can result in fewer costly fermion inversions for a fixed 
acceptance rate. In practice the parameter MSTEPS can be tuned to optimize the update - typically MSTEPS=10. 

Finally, control returns to update () and the updated Hamiltonian Hjiew is computed. A simple Metropolis test 
is used to accept or reject the field configuration at the end of the trajectory. 

4.1. Site, Link and Plaquette type operators 

The bosonic and fermionic fields, and the covariant difference operators living on the hypercubic lattice are asso- 
ciated with various geometric structures such as sites, links and plaquettes. They are implemented in the code using 
various user defined C++ classes: Site_Field, Link_Field, Plaq_Field, Body_Field, etc. They are constructed 
such that they can take values in U{N) or S U{N). They make appearances in the code in many ways and we summarize 
their general structure in the table below: 



Site_Field 


Six) 


^"i^(x) 




Link_Field 








Plaquette_Field 








Body_Field 
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#lnclude "fermion_forces.h" 

void fermion_forces(const Gauge_Field &U, Gauge_Field &f_U, 
const Twist_Fermion &s, const Twist_Fermion &p) 

I 

Lattice.Vector x, e_mu; 
int sites, mu, a, b; 
Umatrix tmp; 
Gauge_Field Udag; 

Udag=Adj(U); 
f_U=Gauge_Field(); 

//contribution to f_U from psi_muDb_mu(U)eta term 
sites=0; 

while(loop_over_lattice(x, sites)) 
{ 

for(mu=0;mu<NUMLINK;mu++) 
{ e_mu=Lattice_Vector(mu); 
tmp=Umatrix(); 

for(a=0;a<NUMGEN;a++) 

I 

for(b=0;b<NUMGEN;b++) 

{tmp=tmp+conjug(p.getS().get(x).get(a))*s.getL().get(x,mu).get(b) 

*Lambda[a]*Lambda[b]*Udag.get(x,mu)-conjug(p.getS().get(x+e_mu).get(a)) 

*BC(x,e_mu)*s.getL().get(x,mu).get(b)*Lambda[b]*Lambda[a]*Udag.get(x,mu);] 

I 

f_U.set(x,mu,f_U.get(x,mu)-0.5*Adj(tmp));) 

I 

sites=0; 

while(loop_overJattice(x, sites)) 
I 

f or(mu=();mu<NUMLINK;mu++) 

{e_mu=Lattice_Vector(mu); 

tmp=Umatrix(); 

for(a=0;a<NUMGEN;a++) 

I 

f or (b=0;b<NUMGEN;b+ +) 

{tmp=tmp+conjug(p.getL().get(x,mu).get(a))*s.getS().get(x+e_mu).get(b) 
*BC(x,e_mu)*Lambda[a]*Lambda[b]*Udag.get(x,mu)- 
conjug(p.getL().get(x,mu).get(a))*s.getS().get(x).get(b) 
*Lambda[b]*Lambda[a]*Udag.get(x,mu);) 

I 

f_U.set(x,mu,f_U.get(x,mu)-0.5*Adj(tmp));) 

I 

sites=0; 

while(loop_overJattice(x, sites)) 

{for(mu=0;mu<NUMLINK;mu++){f_U.set(x,mu,-1.0*Adj(f_U.get(x,mu)));|l 
return; 



Figure 4. A part of the C++ code to compute the fermion force contribution. 
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1 


class Link_Field{ 




2 


private: 




3 


Afield links[SITES][NUMLINK]; 




4 


public: 




5 


Link_Field(void); 




6 


Link_Field(iiit); 




7 


Afield get(coiist Lattice_Vector & 


, const int) const; 


8 


void set(const Lattice.Vector &, 


const int, const Afield &); 


9 


void print(void); 




10 


1; 




11 






12 


Lmk_Field Cjg(const Link_Field &); 




13 


Link_Field operator +(const Link_Field 


&, const Link_Field &); 


14 


Link_Field operator -(const Link_Field 


&, const Link-Field &); 


15 


Linkjield operator *(const double, const Link_Field &); 


16 


Link-Field operator *(const Complex & 


const Link_Field &); 


17 


Complex operator *(const Link_Field & 


, const Link_Field &); 



Figure 5. The Link_Field class with overloading of operators (from utilities .h). 

As an instructive example let us look at the coding details of the Link_Field class. In Fig. 5 we show how the 
Link_Field class is defined along with overloading of basic operators such as '+' and 

We look at the structure of the fermionic term rjD^^^^ on the lattice and the structure of the corresponding fermionic 
operator in the code. On the lattice this fermionic term takes the form 

770^(^^ ^ ^[77(x)^"(A,,(x) + <A/.(x)^i^^77(x)] 

= ^[77°(x)r«©*"VV*(x) + «A^(x)r''©*"V°/,°(x)l , (62) 

where T" are the generators of the gauge group. 

On expanding the lattice covariant difference operators we have 

TTD.iif^ ^ ^[/,"(x)r'(rV^(x)'Z/^(x)-i/;(x-?^)rV'(x-^^)) 

= i[77«(x)(rr*'W^^(x))<A^(x) - Ti\x)[T"a^l{x -7,)T'y^{x -S,)) 

+^^(x)(r''r«'Z/^(x))/7«(x +?^) - iA'(x)(r'"w^(x)r°)77"(x))] . (63) 

In the code we compute the combination Ti(T''^f,(x)T'') as V'^(x)''* and store it as the object Adjoint_Link_Field. 
It is this field that is passed into the functions that require the action of the twisted fermion operator in the inverter. 
Explicitly, the contribution to the operator coming from the term Tr (rjDf^i//^) in the action takes the following form in 
the code: 

+0 . 5*conjug(V. get (x,mu) . get (a,b) ) 

-0 . 5*conjug(V. get (x-e_mu,nm) . get (b,a) ) *BC(x, -e jmu) 

+0 . 5*conjug (V . get (x ,mu) . get (a , b) ) *BC (x , e_mu) 

-0 . 5*conjug (V . get (x ,mu) . get (b , a) ) 

5. Simulation results 

In this section we provide some numerical results obtained through the recent simulations of the two-dimensional 
N ^2 lattice SYM theory t5ni52l . 
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The results we show in this section were obtained using the orbifold prescription for the parametrization of the 
complexified gauge fields ^^{x) on the lattice. The continuum fields ^a{x) are mapped to link fields Ka(ii) living on 
the link between n and in n + /i^, through the mapping: 

U.in) = Jlain) , (64) 

where J?lfl(n) = ^il'i ^'a^' where T' = 1 . . .Nq are the anti-hermitian generators of a U(N) group. Notice though 
that in spite of the appearance of a complex connection the theory only possesses the usual U(N) gauge symmetry. 
Simulations with linear gauge links of this type have been investigated in ifSTl . 

5.1. Eigenvalues of scalars 

The requirement that the lattice theory target the continuum theory as the lattice spacing is sent to zero demands 
vanishing of the fluctuations of all lattice fields and in particular the fluctuations of the trace part of the scalar field 
BJ]. It is also important that the trace mode develops a nonzero expectation value of unity in order that the lattice 
action yield the appropriate kinetic terms in the naive continuum limit. Given the absence of any classical potential 
guaranteeing these features, we find that it is necessary to add a suitable gauge-invariant potential to the lattice theory 
to ensure these conditions holcf*^ In principle, once this mode is regulated one can examine whether this potential 
can be sent to zero in the continuum limit. 

We have added a simple potential term of the following form to regulate the trace mode in the simulation^^ 

This term fixes the vev of the scalar trace mode B° to unity and constrains the fluctuations of the trace mode with 
a quadratic mass term at leading order in the lattice spacing. The remaining traceless fluctuations feel only a soft 
quartic potential. 

SM^^i^Y}SB'if + ... (66) 

X 

Since this U{V) scalar sector decouples in the naive continuum limit this should not break the supersymmetry of the 
remaining S U(N) sector for small enough lattice spacing (indeed all susy breaking terms should vanish as ju — > 0) 



In the C++ code the mass term (65 1 is implemented using 



( 1 . 0/NCDLOR) *Tr (Udag . get (x ,mu) *U . get (x ,mu) ) . real () -1 . 

in action, cpp. The t/(l) mass coefficient fi is denoted by the parameter BMASS and should be held^xeii as we take 
the continuum limit. This implies that the physical mass is taken to infinity in this hmit for any non-zero fi and hence 
that the expectation value of B° is frozen at unity in this limit. 

We rescale all lattice fields by powers of the lattice spacing to make them dimensionless. This leads to an overall 
dimensionless coupling parameter of the form N/(2Aa^), where a - /3/T is the lattice spacing, yS is the physical extent 
of the lattice in the Euclidean time direction and T is the number of lattice sites in the time-direction. The coupling 
A = g^N is the usual 't Hooft parameter. Thus, the lattice coupling 

K = r , (67) 

for the symmetric two-dimensional lattice where the spatial length L - Tp^Note that Aff is the dimensionless physical 
't Hooft coupling measured in units of the area. In these two dimensional simulations, the continuum limit can be 



Notice that our lattice gauge fields are dimensionless and hence contain an implicit factor of the lattice spacing a. 
' ' It was precisely this requirement that led to a truncation of the U (N) symmetry to SU (N) in the original simulations of these theories corTe- 
sponding to a delta function potential for the U(l) part of the field 1 14|. 
'-A potential term of this type was first introduced and tested in [53j. 

"To obtain this dimensionally reduced model from the N = A theory one merely sets the parameters LX = 1, LY = 1 in utilities .h. 
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Q = 4, U(2), BMASS = 1.0, PBC = -1 

0.02 - • t = 0.5 

■ • t=1.0 

• t = 2.0 

0.015 - 



0.01 - 



0.005 - 

! : t I : ; ; : ll 

1 2 3 4 5 6 V 8 9 10 11 12 
Figure 6. Plot showing the average scalar eigenvalue versus the lattice size L in the two-dimensional Q = 4 theory. 

approached by fixing f = and A^, and increasing the number of lattice points L — » oo. We have taken three different 
values for this coupling t = 0.5, 1.0, 2.0 and lattice sizes ranging from L = 2, ■ ■ ■ ,12. In Fig. 6 we show the average 
scalar eigenvalue given by Ul,Ua - I for the Q - 4 model as a function of the lattice size L. This figure confirms that 
as L —> CO we are indeed approaching a continuum limit since the scalar eigenvalues (which contain a factor of a to 
render them dimensionless) are driven to zero. 

5.2. Pfaffian phase/sign problems 

The models we have discussed may encounter an additional difficulty in the context of simulation - the fermionic 
sign problem. After integration over the fermions the effective bosonic action picks up a contribution from the loga- 
rithm of the fermionic Pfaffian Pf(M) which is not necessarily real. Indeed for the supersymmetric lattice constructions 
we described above, M at non zero lattice spacing is a complex operator and one might worry that the resulting Pfaffian 
could exhibit a fluctuating phase e'". Since Monte Carlo simulations must necessarily be performed with a positive 
definite measure the only way to incorporate this phase is through a re-weighting procedure which folds the phase in 
with the observables of the theory. Expectation values of observables derived from such simulations can then suffer 
huge statistical errors which swamp the signal rendering the Monte Carlo techniques effectively useless. 

In Fig. 7 we show results for (| sin(cK)|) as a function of L for the (3 = 4 model with gauge group U{2) (edit 
utilities .h to change number of supercharges). Three values of f = Aff are shown in each plot but the behavior 
is qualitatively similar for all f. We have used the mass parameter controlling the U(l) mode as BMASS = 1. These 
numerical results show that while this model appears to suffer from a sign problem for coarse lattices these effects 
disappear as the lattice is refined and the phase fluctuations are driven to zero as the continuum limit is taken. This is 
consistent with the work reported in 1531 

5.3. Restoration of supersymmetry 

The topological nature of the twisted theory formulated on a torus with periodic boundary conditions can be 
used to show that the partition function of the lattice model is actually independent of the coupling constant. Thus 
derivatives of the partition function with respect to the coupling constant such as the expectation value of the action 
must vanish. Since the fermions enter only quadratically, their contribution can be evaluated simply using a scaling 
argument and thence a simple expression derived for the expectation value of the bosonic action. Thus measurements 
of {S b{'U,'U)) provide us with a check that the scalar supersymmetry has indeed been implemented correctly in our 
codes. Actually, since in practice we use supersymmetry breaking (thermal) boundary conditions (and also employ a 
supersymmetry breaking potential for the scalar t/(l) mode) to do simulations, measuring this quantity provides some 
insight into the magnitude of supersymmetry breaking effects in the theory. 
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Figure 7. Plot showing the average of the sin of the Pfaffian phase a against the lattice size L in the (3 = 4 lattice SYM with gauge group U{2) and 
exponential representation for the gauge links. 



In the case of two-dimensional (3 = 4 theory, we have the expression for the mean action 

<^> = -^ =<^^> = 0, (68) 

OK 

where k is the coupling constant of the twisted action and the last equality follows from the (3-exact nature of the 
twisted theory and shows that the vanishing mean action can be thought of as arising as a consequence of a simple 
S-Ward identity. 

If we integrate out the twisted fermions and the auxiliary field d we find the following expression for the partition 
function of the two-dimensional <3 = 4 theory 



^4NaV/2^-NaV/2 



J IyKD^(e-''^'^'"''"MstkM(^(,^()) , (69) 



where Nq is the number of generators of the gauge group and V is the number of lattice points. The first pre-factor 
arises from the fermion integration and the second derives from the Gaussian integration over the auxiliary field. From 
this we find the following condition on the mean bosonic action as a consequence of the scalar supersymmetry Q: 

(xSb) = ^NcV . (70) 
In Fig. 8 we show the mean bosonic action on the lattice against the lattice size L. The thick solid line represents 



the exact value of the bosonic action given in 70 Clearly the lattice measurements approach the exact result for 
sufficiently small lattice spacing. The deviations that are visible are presumably related to the fact that we have a sign 
problem (these measurements do not incorporate re-weighting) for small L and the simulations are also conducted at 
non zero temperature. We have shown that the sign problem disappears in the continuum limit which is consistent 
with the much better agreement at large L. To recover the true zero temperature result requires in principle that we 
extrapolate our measurements to f — » oo after taking the thermodynamic limit. 



6. Conclusions and outlook 

In this paper we have described in some detail the construction of an object oriented code suitable for the sim- 
ulation of a recently discovered class of lattice field theories possessing exact supersymmetry. The continuum con- 
struction and lattice discretization of (3 = 4, 8, 16 supercharge SYM theories in two, three and four dimensions are all 
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Figure 8. Plot showing the average bosonic action (k(S b + S m)) on the lattice against the lattice size L in the (3 = 4 lattice SYM with gauge group 
C/(2) and exponential parametrization for the gauge links. The thick solid line corresponds to the exact value of the bosonic action. 

covered in detail. The structure of the problem requires the construction of unusual data structures for representing 
the fermions, which is the primary difference between the code described here and more conventional codes suitable 
for simulating QCD. Nevertheless the basic algorithms employed (RHMC and multi-mass CG solvers) are borrowed 
directly from lattice QCD and adapted to the problem at hand. We verify the correctness of the resultant code by 
showing results from simulations of the two-dimensional SYM model. Acceleration of this code can be achieved by 
off-loading the linear solver calculation to a GPU card - we refer the reader to ||54| for details. It is also possible 
to parallelize the code with suitable distributed libraries layered over MPI f55l and work in both these directions is 
ongoing. 
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Appendix A. Installation of the program 

It is very easy to perform the installation and execution of SUSYXATTICE. Below we provide the necessary 
steps on Unix or Linux systems. 

1 . Download the code from CPC Program Library and unpack it. 

2. Change the directory to SUSYXATTICE. 

3. Compile the code (g++ -O *.cpp -o SUSYXATTICE -llapack -Iblas). 

4. Modify the input parameters located in file parameters 

5. Type ./SUSYXATTICE >& log & to run the code. 

The authors have tested the code on Linux machines. After slight modifications of above steps the code may be 
installed on other machines. 

The output of the code produces the following files in the running directory: 

1. cgs : Average number of conjugate gradient (CG) iterations. (See MCG_solver . cpp). 
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2. conf ig: File to read in containing the site, link and plaquette field configurations from a previous run. (See 
read_in. cpp.) 

3. corrlines: Correlation function between temporal Polyakov lines as function of spatial separation (See 
corrlines . cpp.) 

4. data: Boson (1st column) and fermion (2nd column) contributions to the total action. (See measure . cpp.) 

5. dump : Site, link and plaquette field configurations stored as ASCII (See write_out . cpp.) 

6. eigenvalues : Eigenvalues of IC^aggeraixyilaix) N real numbers for each lattice point x and direction a (See 
measure . cpp.) 

7. hmc_test: e"'"'''"'^ from HMC test (See update . cpp.) 

8. lines_s : Spatial Polyakov line. (See measure . cpp.) 

9. lines_t : Temporal Polyakov line. (See measure . cpp.) 

10. log: Log file. 

11. loops: Wilson loops. (See loop . cpp.) 

12. scalars: Tr(f/' f/) (See measure . cpp.) 

13. ulines_s : The spatial Polyakov line computed using the unitary part of the link (See measure . cpp.) 

14. ulines_t : The temporal Polyakov line computed using the unitary part of the link (See measure . cpp.) 

Appendix B. The list of files in SUSY_LATTICE library 

We list the files included in SUSYXATTICE library with a brief description of their purpose. 

1. action, cpp: Compute the total action - fermionic and bosonic. 

2. corrlines . cpp : Finds the traced product of the link matrices at various lattice sites. 

3. evolve_f ields . cpp : Leapfrog evolution algorithm. Also stores the fermion and boson forces for the next 
iteration. 

4. f ermion_f orces . cpp : Computes the fermion kick to gauge link force. 

5. force, cpp: Bosonic and pseudo-fermionic contribution to the force. 

6. kinetic_energy . cpp : Computes the kinetic energy term in the Hamiltonian. 

7. line, cpp: Computes the Polyakov lines. 

8. loop . cpp : Computes the Wilson loops. 

9. matrix . cpp : Builds the fermion matrix (sparse and fuU forms) and also computes the Pfaffian of the fermion 
operator 

10. MCG_solver . cpp : multi-mass CG solver needed for RHMC alg. 

11. measure . cpp : Performs measurements on field configurations. Writes out scalar eigenvalues, Polyakov/Wil- 
son loops and the action. 

12. my_gen.cpp: Computes 5 t/(A^) generator matrices. 
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13. obs . cpp : Computes fermion and gauge actions. Also returns the unitary piece of the complex link field. 

14. read_in . cpp : Reads in the previously generated field configurations - file conf ig 

15. read_param . cpp : Reads in the simulation parameters from a data file called parameters 

16. setup . cpp : Contains the partial fraction coefficients necessary to represent fractional power of fermion oper- 
ator - used by Remez algorithm. 

17. sym . cpp : The main program - performs warm up on field configurations and commences measurement sweeps 
once the configurations are warmed up. 

18. unit . cpp : Extracts the unitary piece of the complex gauge links. 

19. update. cpp: Updates the field configurations based on HMC test. 

20. utilities . cpp : Utility functions. Contains constructors for site, link, plaquette fields, gauge fields, twist 
fermions etc. Edit to change number of supercharges and size of lattice dimensions. 

21. write_out . cpp : Writes out the values of gauge and twist fermion fields on to a file called dump. 

Appendix C. A sample input parameter file for SUSY_LATTICE 

This is a sample input parameter file called parameters located in the SUSY_LATTICE folder. 
1008(9 50 10 0.5 1.0 0.02 0.0 

SWEEPS THERM GAP LAMBDA BETA DT ALPHA READIN 
There are the definitions of the parameters: 

1 SWEEPS : Total number of Monte Carlo time steps intended for taking measurement steps. 

2. THERM : Total number of Monte Carlo time steps intended for thermalizing the field configurations. 

3. GAP: The gap between measurement steps. 

4. LAMBDA : The 't Hooft coupling. 

5. BETA : Inverse temperature. 

6. DT : The time step put in the integrator for leapfrog evolution. 

7. ALPHA : A supersymmetric mass (deformation) parameter. 

8. READIN : Determines whether to read in the previously generated field configurations or not. The program will 
read in the previous configurations if READIN is set to 1. 
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PROGRAM SUMMARY 

Manuscript Title: An object oriented code for simulating supersymmetric Yang-Mills theories 

Authors: Simon Catterall and Anosh Joseph 

Program Title: SUSY_LATTICE 

Journal Reference: 

Catalogue identifier: 

Licensing provisions: None 

Programming language: C++ 

Operating system: Any, tested on Linux machines 

Keywords: Lattice Gauge Theory Supersymmetric Yang-Mills Rational Hybrid Monte Carlo Object Oriented Programming 
PACS: 11.15.Ha, 12.60.Jv, 12.10.-g, 12.15.-y, 87.55.kd, 87.55.kh 
Classification: 1 1.6 Phenomenological and Empirical Models and Theories 
Nature of problem: 

To compute some of the observables of supersymmetric Yang-Mills theories such as supersymmetric action, Polyakov/Wilson 
loops, scalar eigenvalues and Pfaffian phases. 
Solution method: 

We use the Rational Hybrid Monte Carlo algorithm followed by a Leapfrog evolution and a Metroplois test. The input parameters 

of the model are read in from a parameter file. 

Restrictions: 

This code applies only to supersymmetric gauge theories with extended supersymmetry, which undergo the process of maximal 
twisting. (See Section|2]of the manuscript for details.) 
Unusual features: 
Running time: 

From a few minutes to several hours depending on the amount of statistics needed. 
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